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ABSTRACT 

A method of calculation is presented that allows the 
simulation of the time-dependent three-dimensional mo- 
tion of thin liquid layers on solid substrates for systems 
with finite equilibrium contact angles. The contact an- 
gle is a prescribed function of position on the substrate. 
Similar mathematical models are constructed for sub- 
strates with a pattern of roughness. Evolution equations 
are given, using the lubrication approximation, that in- 
clude viscous, capillary and disjoining forces. Motion 
to and from dry substrate regions is made possible by 
use of a thin energetically- stable wetting layer. We sim- 
ulate motion on heterogeneous substrates with periodic 
arrays of high contact-angle patches. Two different prob- 
lems are treated for heterogenous substrates. The first 
is spontaneous motion driven only by wetting forces. 
If the contact-angle difference is sufficiently high, the 
droplet can find several different stable positions, de- 
pending on the previous history of the motion. A second 
simulation treats a forced cyclical motion. Energy dissi- 
pation per cycle for a heterogeneous substrate is found 
to be larger than for a uniform substrate with the same 
total energy. The Landau-Levich solution for plate re- 
moval from a liquid bath is extended to account for a 
pattern of roughness on the plate. 


1. INTRODUCTION 

Wetting and capillary considerations during the slow 
motion of liquids on solid substrates are important in 
both the technological and natural worlds. Applications 
include the spreading behavior of liquid coatings, as 
well as flows in oil reservoirs, chemical reactors and 
heat exchangers. Some biological application areas are 
motions in the tear film on the cornea of the eye, flows 
on liquid covered membranes in the lungs, and the gen- 
eral area of cell motility, where cell motions have many 
of the features of inert liquid droplets [ 1 ]. The spreading 
properties of agrochemicals such as pesticides and in- 
secticides are important determinants of their effective- 
ness. Capillary forces are often dominant in problems 
of small physical dimensions; however they can also be 
important for large-scale phenomena when body forces 


are very weak, as in the microgravity environment of 
orbiting satellites in space. 

Application of a wettability pattern or a pattern of 
roughness to the walls of vessels containing liquid is a 
possible strategy for controlling the position of the liq- 
uid in the microgravity environment of orbiting space 
vehicles. In the absence of gravity, liquid in a partially- 
filled container can assume many different configura- 
tions. Small imposed accelerations can cause large dis- 
placements of the liquid which is undesirable for a va- 
riety of reasons. Contact line dynamics will play an 
important role, either by “pinning” the location of the 
liquid on the wall or by dissipating the kinetic energy 
imparted to the liquid. The present results suggest that 
different wall wettability patterns may accomplish one 
or the other of these objectives. 

We will briefly outline mathematical and numerical 
procedures for liquid motions on chemically heteroge- 
neous or roughened substrates. These employ the long- 
wave or “lubrication” approximation. More complete 
information can be found elsewhere [2,3,4]. 


2. THE MATHEMATICAL MODEL 
Integral mass conservation is 

h( = -V • Q + Wi(x,y,t) . (2.1) 

Here h is the thickness of the liquid layer and Q is the 
two-dimensional areal flux vector defined as 



where u and v are velocity components in the x and y di- 
rections respectively. Likewise V is a two-dimensional 
operator with respect to x and y. w; is a local injection 
rate that is an input function of substrate position and 
time. 

Under the assumptions of lubrication theory [5,6], 
i.e. the motion is sufficiently slow that inertial forces 
may be neglected and the free surface is inclined at a 
small angle relative to the substrate, the momentum equa- 
tion for an assumed Newtonian liquid may be integrated 
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to yield 

Q = -±{h*Vp) . (2.2) 

Here // is the viscosity and the pressure p is independent 
of the normal coordinate z. The no-slip condition for the 
velocity has been applied on the substrate and the liquid 
free surface is stress- free. The pressure in the liquid is 
given by 


p=-< 7 /v-n n (2.3) 


and the pressure above the liquid is taken to be zero 
without loss of generality. The first term on the right 
is an approximation to the free-surface curvature when 
the surface slope is small. The error in this curvature ap- 
proximation is proportional to the square of the surface 
inclination, a is surface tension and the so-called dis- 
joining pressure is given by the two-term model [ 2 , 3 * 7 ] 



D and the exponents n and m are positive constants with 
77 > m > 1. The local disjoining energy density 


e {d \h) = - 



(2.5) 


has a single stable energy minimum at the thin wetting- 
layer thickness h = /i*. A local force balance near an 
apparent contact line gives 

acosO c = a — c^(oo) ( 2 - 6 ) 


which is the disjoining- model equivalent of the Young- 
Laplace equation. 

The equilibrium contact angle is taken to be an ar- 
bitrarily prescribed function of position and the evolu- 
tion equation for flow over a heterogeneous substrate 
becomes 


Each pressure component on the right of equation 
( 2 . 3 ) may be identified with an integrated energy com- 
ponent. The free-surface energy is proportional to the 
area of the liquid surface and is given by 

(2.8) 

where A is the total area of the substrate, and 7 is the 
angle between the normal to the surface and the normal 
to the substrate. The total disjoining energy is 

E w = JJ ? ld) dA. (2.9) 

The global energy change equation is 

£<") _ _(£(*) + E (d) ) + ir (2.10) 

which is the statement that the rate of viscous working 
is equal to the rale of decrease of the stored, or poten- 
tial, energy components plus the rate of working on the 
system. Here 

I v = JJpwjdA (2.11) 

may be recognized as the rate of working on the system 
by means of injection w t . For cyclical motions, the drop 
shape returns to its original configuration and there is no 
change, over a cycle, in or In that case the 
total viscous work done is the time integral of IF or 

j £ p dV 

where F is drop volume, rfF = uu dA and the special 
integral sign denotes a full cycle. 

For the quasi-three dimensional problems treated here, 
numerical solutions use finite difference methods and an 
alternating direction implicit (ADI) technique is imple- 
mented. Developed originally for second-order elliptic 
and parabolic systems [ 8 ], ADI uses alternating sweeps 
in each direction and only a banded system of equations 
needs to be solved to update the discrete set of hj j val- 
ues. 


) + w> 

(2.7) 

It is possible to augmented this equation with gravity 
and other force terms as additional driving mechanisms. 

[3] 


3. VALIDATION 

The model has been calibrated by comparison with 
several experiments. Axisymmetric spreading of a droplet 
on a high-energy [i.e. small contact angle] uniform sub- 
strate has been measured by several investigators and 
general agreement with Tanner’s Law [9] 
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has been verifcd. Here A is drop area, V is volume and 
t is time. It can be shown that the solution for h(r. t) is 
self-similar and the drop central height must also follow 
a power law in time [3]. Figure 1 shows a compari- 
son of the similarity solution assuming a wetting layer 
of 30 nm thickness which agrees with experimental re- 
sults [10]. Also shown is an axisymmetric simulation 
on a uniform low-energy substrate with a finite contact 
angle. It is seen that the power- law is followed until 
the drop has almost stabilized at a dimensionless cen- 
tral height h c » 1. This spreading behavior, with fi- 
nite contact angle, has been observed experimentally by 
Zosel [11]. Another simulation, shown as symbols in 
the Figure, assumes the substrate has a pattern of wet- 
tability as in Fig. 5. The substrate area- average wetta- 
bility is the same as for the low-energy uniform case. 
It can be seen that the spreading is unaffected by the 
contamination pattern until the droplet slows and ulti- 
mately stops spreading. Three-dimensional simulations 
are limited to thick precursor layers because of the need 
to adequately resolve details in the contact region. Too 
thick a precursor layer results in a relatively weak speed 
discrepancy; this discrepancy is in accordance with the 
predicted inverse logarithmic dependence on the precur- 
sor thickness [12]. 



Figure 1 : Drop central height h c versus time, showing 
the similarity range for axisymmetric spreading. The 
initial drop height is 10. The straight line on this log- 
log plot is h c = 0.75* 1 / 5 which matches experiments. 
The solid curve is the calculated behavior for finite con- 
tact angle on a uniform substrate, while the lines-points 
curve is calculated for the square wettability pattern. 

We have performed a laboratory experiment to com- 
pare with drop spreading simulations [3]. A 26 fA drop 
of glycerin was placed near the center of a cross of 1 



Figure 2: Video images taken from an experiment where 
a glycerin drop is placed near the center of a cross of 0. 1 
cm Teflon tape on a glass slide. These pictures may be 
compared with the simulation shown in Fig. 3. 



Figure 3: Contour plots from a computer simulation of 
break-up of a liquid drop placed on a “cross” of high- 
contact- angle material. 
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Figure 4: Energy components versus time from the 
drop-on-cross simulation. 


mm Teflon tape that had been fixed to a horizontal glass 
slide. Wetting forces cause the drop to break up into un- 
equal fragments as shown in Fig. 2. Simulation results, 
in Fig. 3, show detailed agreement with the experiment. 
However time-scale corrections need to be applied be- 
cause (i) the simulation precursor layer is overly large, 
and (ii) the contact angles in the experiment are beyond 
the range of quantitative validity of the small-slope lu- 
brication approximation. Energy component variation, 
from the simulation, is given in Fig. 4. The Bond num- 
ber Do = pgR 2 /a = 2.5 using the stabilized radius R 
for the drop on glass. The gravitational effect is quite 
minor, however, as can be seen in the Figure. Note that 
the motion proceeds in a “jerky” manner that is charac- 
teristic of capillary driven motions on nonuniform sub- 
strates. Each drop disconnection is reflected as a rapid 
decrease in the substrate energy. 


4. HYSTERETIC MOTIONS ON HETEROGENEOUS 
SURFACES 

Simulations have been performed to investigate the 
dissipation of energy when a liquid moves on mixed- 
wettable substrates. Again we consider the motion of a 
drop, both in spontaneous motion, driven only by wet- 
ting forces, and also in a periodic forced motion. Sev- 
eral periodic patterns of wettability have been consid- 


(a) 



(b) 



Figure 5: The substrate wettability patterns used in this 
study. The dark regions have larger values of contact an- 
gle than the surrounding field, (a) Isolated spot pattern; 
(b) Checkerboard pattern. 


cred, two of which are shown in Fig. 5. The equilibrium 
contact angle 0 e is larger on the dark patches than on the 
surrounding field. Numerical procedures for generating 
these and other patterns, including slight smoothing at 
patch boundaries to maintain derivative continuity, are 
given in [2]. 

Figure 6 shows two different stable drop shapes ob- 
tains by allowing a drop to advance outward (left) and 
recede (right) from a starting profile that was cither steeper 
or shallower than the equilibrium shapes. The contact 
angle on the isolated patches is about four times larger 
than on the field. Effective advancing and receding con- 
tact angles may be calculated using averages of local 
values; the effective contact angle ratio for this case is 
e a jO r - 1.14. This ratio is a measure of contact an- 
gle hysteresis. Figure 7 is an instantaneous picture dur- 
ing the receding simulation. It shows transient interior 
dewetting due to an instability near the receding contact 
line. 

A schematic diagram for a notional experiment to 
explore energy dissipation in forced motions in shown 
in Fig. 8. By pumping liquid in and out via a syringe, 
the droplet can be forced to periodically traverse the 
wettability pattern. Pressure-volume plots from numer- 
ical simulation are shown in Fig. 9 where motions on 
a checkerboard wettability pattern are compared with 
a uniform substrate of the same average energy. The 
area of each hysteresis loop is the input work required 
to drive a cycle of the motion. Additional dissipation 
results from the presence of the pattern, the effect be- 
ing more important at low speeds or long cycle times. 
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Figure 6: Contour plots of final static shapes for ad- 
vanced (left) and receded (right) drops in spontaneous 
motion. The wettability pattern is also shown. The dif- 
ference between the two cases is a measure of contact 
angle hysteresis. 



Figure 7: A frame during recede. A dc- wetting insta- 
bility leads to interior dry patches near the receding pe- 
riphery. This is a transient effect and the drop continues 
to recede until the final configuration, shown in Fig. 6, 
is attained. 



Figure 8: Schematic diagram of a notional experiment 
to investigate energy dissipation in a cyclic motion of a 
drop. 



Figure 9: Dimensionless pressure P versus volume V 
showing hysteresis loops in cyclical motion. The solid 
line is for a uniform substrate, dashed line is for a 
checkerboard wettability pattern. The symbol curve is 
for the non-dissipative quasi-static theory with equation 
P = 27r I / 3 r~ 1 / 3 . 

For comparison, a quasi-static non-dissipative curve, for 
motion on a uniform substrate, is also shown. 


5. MOTIONS ON A SURFACE WITH A PATTERN 
OF ROUGHNESS 

It is possible to include a roughened substrate in the 
lubrication formulation. It may be shown [4] that the 
evolution equation needs only be modified by inclusion 
of the substrate shape function in the “permeability,” u 
e. the factor of proportionality between the flux and 
the pressure gradient. It is often convenient to allow 
the substrate to move with time while the computational 
window is fixed to the liquid free surface. 

A model problem is an extension of the well-known 
Landau-Levich [13] problem for the withdrawal of a 
moving plate from a bath of liquid. We consider the 
plate to have a periodic pattern of roughness or “cells” 
as shown in Fig. 10. The plate moves with constant 
speed V to the right while the liquid meniscus is pinned 
at the left end of the computational window. While 
the Landau-Levich result is time independent, here the 
problem becomes time periodic with a period equal to 
the cell passage time. All dimensional constants can be 
absorbed by scaling. We let ho be a measure of coat- 
ing thickness, such as the cell depth and substrate coor- 
dinates (x, y) are made dimensionless using the length 
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Figure 10: A periodic square array of cells on the sub- 
stratc. 



Figure 1 1 : A liquid surface during meniscus withdrawal 
for the cell pattern in Fig. 10. Note how the meniscus 
’scrapes out’ the cells. 



The time unit is 

r* = (5.16) 

u 

The long-wave model equation is 

U ~ - v • (sVV 2 /,) (5.2) 

at ox 


A function Jn(x,yJ) represents the moving substrate 
and the permeability s in (5.2) is simply s = (h — h \ ) 3 . 

As seen in Fig. 1 1 , the downward pressure of the 

meniscus acts so as to “scrape” liquid from the cells. 
A certain residual fraction remains in each cell, deter- 
mined primarily by the cell geometry. Much greater 
work is required to move the plate compared to a smooth- 
wall case. The problem is applicable to damping of peri- 
odic motions where the wall roughness is a surrogate for 
the chemical heterogeneity discussed above. Further de- 
tails, including terrestrial applications in the coating and 
printing industry, where engraved or “gravure” rollers 
are used, are given in [4]. 
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